February 2, 2008 7:0 WSPC/INSTRUCTION FILE ee 



International Journal of Theoretical and Applied Finance 
© World Scientific Publishing Company 



THE EXPONENT EXPANSION: AN EFFECTIVE 
APPROXIMATION OF TRANSITION PROBABILITIES OF 
DIFFUSION PROCESSES AND PRICING KERNELS OF 
FINANCIAL DERIVATIVES 



LUCA CAPRIOTTI 

Global Modelling and Analytics Group, Credit Suisse 
One Cabot Square, London, EH 4QJ> United Kingdom 
luca. capriotti@credit-suisse. com 

Received (19 September 2005) 
Revised (18 January 2006) 

A computational technique borrowed from the physical sciences is introduced to obtain 
accurate closed-form approximations for the transition probability of arbitrary diffusion 
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1. Introduction 

Continuous-time diffusion processes are at the basis of much of the modeling work 
performed every day in Finance and Economics. Indeed, since Bachelier's pioneering 
work on the application of probability theory to the dynamics of stock prices [1], 
many economic variables subject to unpredictable fluctuations, have been modeled 
by stochastic differential equations of the form 

dY t = Vy{Y t )dt + <Jy(Y t )dW t . (1.1) 

Here ^ y {y) is the drift, describing a deterministic trend, and <J v (y) > is the 
volatility function, describing the level of randomness introduced by the Wiener 
process (i.e., white noise), dWt- The main reason for the popularity of this class of 
models is probably that in continuous time one can perform analytic calculations 
using the instruments of stochastic calculus, and the powerful framework of partial 
differential equations. In particular, for the few cases for which the process (1.1) is 
exactly solvable, one can derive closed-form solutions for the associated transition 
probability. The latter contains all the statistical properties of the financial quantity 
modeled by the diffusion, and can be exploited in a variety of ways, including 
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the derivation of no-arbitrage prices for financial derivatives in complete markets. 
The milestone results derived by Black, Scholes and Morton [2], Cox, Ingersoll and 
Ross [3] , or Vasicek [4] , are among the most significant examples of the amount of 
progress in Economics that has been done using integrable continuous-time diffusion 
processes. 

Nonetheless, an accurate description of the market observables requires in gen- 
eral more sophisticated models than those for which an analytic solution is available. 
These are usually tackled by means of numerical schemes ultimately relying on a 
discretization of the diffusion, obtained by replacing the infinitesimal time dt with 
a finite time step, At. The approximate results obtained in this way become exact 
only approaching the limit At — ► 0, and this can be done with some computational 
effort. 

In this paper, wc utilize the exponent expansion - a technique introduced in 
chemical physics by Makri and Miller [5] - to derive a short-time approximation 
of the transition probability of the diffusion process (1.1). The aim is to obtain 
an analytic approximation which is as accurate as possible for a time step At as 
large as possible. On one hand, this allows one to derive approximations of financial 
quantities that are very accurate even for sizable values of the time step. On the 
other, it allows a reduction of the computational burden of numerical schemes as 
the limit At — > can be achieved with larger time steps, i.e., with less calculations. 

The possibility to use Makri and Miller's technique to derive approximations of 
the transition probability was originally hinted by Bennati et al. in Ref. [6]. Here 
we explore this possibility, giving derivations for a generic diffusion process with 
state-dependent drift and volatility, and we study the reliability of the exponent 
expansion by applying it to several diffusion processes of financial interest. 

Through the exponent expansion, the transition probability is obtained as a 
power series in At which becomes asymptotically exact if an increasing number of 
terms is included, and provides remarkably accurate results even truncating it to 
the first few (say, n = 3) terms. Two derivations are offered, the first by means of 
Kolmogorov's forward equation [7] (Sec. 2), and the second introducing a slightly 
different formalism (Sec. 2.1). The latter, once the problem is formulated in terms of 
Feynman's path integrals [8] , allows the generalization of the exponent expansion to 
the calculation of the pricing kernel of financial derivatives whose underlying follows 
the considered diffusion. This allows in turn the derivation of simple approximations 
for the price of such contingent claims (Sec. 3). In Sections 2.2 and 3.1.2, we illustrate 
the exponent expansion through the application to the Vasicek, the Cox-Ingersoll- 
Ross, and the Constant Elasticity of Variance models, and in Section 4 we discuss its 
application to Monte Carlo calculations within the path integral framework [9,10,11, 
12]. Finally, we draw our conclusions and we discuss further possible developments 
in Section 5. 
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2. Exponent expansion and transition probabilities of diffusion 
processes 

In this Section, we derive the exponent expansion for the transition probability, 
p v (y, At\y n ), giving the likelihood that the random walker following the process 
(1.1) ends up in the position y at time t — At, given that it was in y at time 
t = 0. Although the derivation can be generalized to the case where drift and 
volatility have an explicit time dependence, here for simplicity we will restrict the 
discussion to the time homogeneous case. In order to make the derivation easier, 
it is convenient to transform the original process in an auxiliary one, say X t , with 
constant volatility a x . Following Ait-Sahalia [13], this can be achieved in general 
through the following integral transformations 

* dz 



X t = 7 (y t ) = ±a x 



<7y{z) 



(2-1) 



where the choice of the sign is just a matter of convenience depending on the specific 
problem considered. The latter relation defines a one to one mapping between the 
x and y processes as the condition <J y {z) > ensures that the function x = 7(2/) 
defined by (2.1) is monotonic, and therefore invertible. a A straightforward appli- 
cation of Ito's Lemma [7] allows one to write the diffusion process followed by X t 
as 



with 



Px{x) 



dX t = p x {X t )dt + a x dW t 



^(7- H*)) 2 dy 



-^(7" 1 W) 



(2.2) 



(2.3) 



where y = 7 _1 (x) is the inverse of the transformation (2.1). Finally, the transition 
probability for the x and y processes are simply related by the Jacobian associated 
with (2.1) giving 

Px{l{y),At\x ) 



p v (y,At\y ) 



(2.4) 



In order to find an expression for the transition probability associated with 
Eq. (2.2) which is accurate for a time At as long as possible, we make the following 

ansatz: 



p x (x, At\x ) 



1 



exp 



(x - x ) 2 
2alAt 



W(x,x ,At) 



(2.5) 



Such transition probability must satisfy the Kolmogorov forward (or Fokkcr- 
Planck) equation [7]: 



d t p(x,At\x ) 



1 



-d x p x (x) + ^ 2 x df 



p x (x,At\x ) 



(2.6) 



For a discussion of the regularity conditions on the drift and volatility functions see e.g., Ref. [13]. 
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This implies in turn that the function W(x,xo,i) follows the equation: 



1 



1 



X — Xq 



t W = -^ x d x W+-a 2 x diW--ai {d x W)* + d xllx — A/ 
Expanding the function W(x,xo,t) in powers of At 7 

oo 

W(x,x ,At) = J2 w n(x,x )At n , 

n=0 



d x W+!±) . (2.7) 



(2.8) 



substituting it in Eq. (2.7), and equating equal powers of At leads in a straightfor- 
ward way to a decoupled equation for the order zero in At giving 



1 r 

W (x,x ) = g / dzfi x (z) , 

and to the following set of recursive differential equations: 

~1 



(2.9) 



2 a A - Vxdx 



(n + l)W n+1 = -(x - x )d x W n+1 + 

_^ m—n 



(2.10) 



m=0 



In particular, for n = 0, 1, 2 Eqs. (2.10) read: 
Wi{x,xq) = -{x - x Q )d x W 1 (x,x ) + 



2^2 M 35 ) 2 + ^d x fi x (x) 



2W 2 (x,x ) = -(x - x ) d x W 2 (x,x ) + ^a 2 x d x Wi{x,x ) , 
3W / 3 (x,a;o) = -(x - x ) d x W 3 (x, x ) + ^ald x W 2 (x, x ) , 



(2.11) 
(2.12) 

(2.13) 



The differential equations above (2.10) are all first order, linear and inhomogeneous 
of the form 



nW n (x,x ) = -{x - x )d x W n {x,x ) + A n -i(x,x Q ) , 



(2.14) 



where A„_i(i,io) is a function that is completely determined by the first n — 1 
relations. It can be readily verified by substitution and integration by parts that 
the solution of (2.14) reads 

W n (x,x )= [ d^fX-ifo + (*-*(>)£, *o) . (2.15) 
Jo 
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This, for n = 1,2,3, after some manipulations, gives: 

Wifoaro) = 4- f dzV cS (z) , (2.16) 



Ax _ 



Xo 



W 2 (x,x ) = ^^[V eS (x)+V eS (x )-2W 1 (x,xo)] , (2.17) 



W 3 (x,x ) 



2Ax 4 
3a 2 



Ax f dzV eS (z) 2 -U dzV eS (z) 



^ 2 W 2 {x,x ) + ^^[d x V eS {x)-d x V eS (xo)] , (2.18) 

where Ax = x — xo, and, for reasons that will be clearer in the next Section, we have 
also introduced the 'effective potential' as the following quantity with dimension 

time^ 1 : 

V eS (x) = -^^ x (xf + ^d x fi x (x) . (2.19) 

At this time, we just note that the first order correction can be rewritten as 

1 f At 

W 1 (x,x ) = —J dtV cS (x + tAx/At) , (2.20) 

leading to the interpretation of this term as a time-average of the effective potential 
over the straight line, constant velocity (Ax/ At) trajectory between xq and x. 
Similarly, the leading term in Ws(x,xq) (i.e., the one proportional to the lowest 
power of the volatility) is proportional to the variance of the effective potential over 
the same trajectory. Finally, we observe that the corrections W n (x, xq) are well 
defined in the limit Ax — > 0. In particular, for n = 1, 2, 3 it is not difficult to show 
that 

lim Wi{x,x ) = V eS {x ) , (2.21) 



*xo 



2 



lim W 2 (x,x ) = 8 2 x V e n{x) , (2.22) 

x^x a 12 

lim W 3 (x,x ) = -§ (dMx)) 2 + £ d 4 x V cS (x) . (2.23) 
x^x„ 24 240 

The form of the trial transition probability represents the main difference of 
the present approach to the one proposed in Ref. [13], which is otherwise very 
similar in spirit. In fact, the latter expands in powers of At the exponential 
exp [— W(x, xq, At)] rather then just the exponent, as we do here, instead. As it 
will be shown explicitly in the following, the present choice gives rise to a distinct 
approximation scheme for n > providing generally a similar level of accuracy 
but remarkably simpler mathematical expressions. This is because, by keeping the 
exponential form of the ansatz, one formulates a guess which is closer to the exact 
one. The latter, can be expected to have an exponential form in order to satisfy 
the Chapman-Kolmogorov property of Markov processes [7] . In addition, the expo- 
nential choice of the ansatz automatically enforces the positive definiteness of the 
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transition density which is not granted in the approach of Ref. [13]. In fact, as it 
will be shown explicitly in Sec. 2.2, in the limit of very small volatility (a x — > 0), 
when the effect of the noise disappears and transition density converges to a Dirac's 
5 distribution, the expansion of Ref. [13] breaks down as the transition probability 
becomes negative. On the contrary, the exponent expansion remains well defined 
and accurate also in this limit. Indeed, the first terms of the expansion in At can 
be also derived through a small volatility expansion of the transition density, as it 
will be discussed in Sec. 3.1. 

Similarly to the expansion developed in Ref. [13] the exponent expansion has in 
general a finite convergence radius which is a decreasing function of the volatility. 
As it will be shown in the following, for the values of volatilities and At relevant for 
financial applications the exponent expansion turns out to be very accurate even 
when truncated to the first few terms. 



2.1. Alternative derivation 



The term Wq(x, xq) in the exponent expansion is somewhat different from the higher 
order terms. In fact, it is defined by Eq. (2.9) which is decoupled from the recursive 
system (2.10). Indeed, it is possible to obtain the same result for the exponent 
expansion by expressing the transition density as 



Px (x,At\x ) = e- w °( x < x °^ Px (x,At\x ) 



(2.24) 



and looking for an approximate expansion of the form (2.5) for <& P:c (x, At\xo). It is 
easy to show by direct substitution in the forward Kolmogorov equation (2.6) that 
$ Px (x, At|a;o) is the solution of 

H x <^ Pl {x, At\x ) = -d t <5> Px (x, At\x ) , (2.25) 

where Ti x is the "Hamiltonian" differential operator 

H x = -^-d 2 x + V cS {x) , (2.26) 

and V e g(x) is the effective potential of Eq. (2.19). As a result one can equivalcntly 
derive the exponent expansion by substituting in Eq. (2.25) the following trial func- 
tion 



$ Px (x, At\x ) = 



^h^At 



exp 



[x - Xq)' 
2a?. At 



-J2w n (x,x )At n 



71=1 



(2.27) 



which does not contain the term Wq{x, xq). This observation will be used in Section 
3 to generalize the exponent expansion to the pricing kernel of financial derivatives. 



2.2. Applications 

The application of the exponent expansion to a generic diffusion process of the form 
(1.1) is rather straightforward and reduces to the calculation of one dimensional 
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Fig. 1. Accuracy of the exponent expansion for the transition density of the Vasicek model (2.28). 
for a = 0.0717, b = 0.261 and xq = 0.1. Panel (a): comparison between the exponent expansion 
(stars) and the approximation of Rcf. [13] (squares) for a = 0.02237 and At = 0.5. At order 
zero the two schemes are identical. The uniform error of the Euler approximation is also reported 
for comparison. Panel (b): maximum relative error for a = 0.02237 as a function of At: n = 
(continuous), n = 1 (dotted), n = 2 (long dashed) and n = 3 (short dashed). The inset is an 
enlargement of the 5-10 years region. Panel (c): comparison between the exponent expansion 
(stars) and the approximation of Ref. [13] (squares) in the regime of low volatility (<r = 0.01), for 
At = 0.5. 



integrals. In this Section, we illustrate this procedure for a few test cases, namely for 
the Vasicek [4] , the Cox-Ingersoll-Ross [3] , and the Constant Elasticity of Variance 
[14] diffusion processes. We will compare the results of the exponent expansion with 
the exact results available in literature, and with the approach of Ref. [13]. 



2.2.1. Vasicek diffusion 

We first consider the Ornstein-Uhlenbeck diffusion proposed by Vasicek [4] as a 
model for the short-term interest rate: 



dX t = a(b - X t )dt + adW t , 



(2.28) 
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where a, b, and a are positive constants representing the mean-reversion level, the 
velocity to mean reversion, and the volatility, respectively. This model is integrable 
and the corresponding probability density function is Gaussian: 



p cx (x,At\x ) = 



1 



(2^(T 2 ) 1 /2 



exp 



[(.x - a)e- aA * - (x - a)] 



121 



2a 2 



with 



^ _ g -2aAt 



2a 



(2.29) 



(2.30) 



The exponent expansion of the Vasicek model can be easily derived using 
Eqs. (2.9), and (2.16-2.18) with the effective potential, Eq. (2.19), 

a 2 (b-x) 2 a 
Wx) = 2 ' 

and gives, 



(2.31) 



p x (x, At\x ) = 



1 



V2na 2 At 



exp [- 



(x - x ) 5 



W (x,x ) 



2a 2 At 

- Wi{x,x )At - W 2 (x,x )At 2 - W 3 (x,x )At 3 ] 
up to the third order in At (n = 3). Here 

a(x - b) 2 - a(xo - b) 2 



W (x, x ) 
Wi(x,x ) 
W 2 (x,x ) 
W 3 (x,x ) 



a 

6a 2 



((x 



2Ax 2 1 



2Ax 3 



2o 2 

- bf + {xo - b) 2 + (x - b)(x - b)) 
V cS (x) + V cS (x ) - 2W 1 (x, x )] , 
^- 4 [( x -bf-(x -bf] 



a 
2 



(2.32) 

(2.33) 
(2.34) 
(2.35) 



a -Ax-^- 2 [(x-bf-{x»-bY 

r 2 



+ 



4 6<7 2 



: {W l (x,x )Y-—^W 2 {x,x {) ) + 



2 2 



(2.36) 



2Ax 2K iy ' v " Aa J ir ' 4Ax 2 
with Ax = x — xq. It is interesting to note that the approximate transition proba- 
bility obtained with the present approach reproduces exactly the expansion of the 
exact transition density Eq. (2.29) at the same order. 

The fast convergence of the approximation scheme is illustrated in Fig. 1. Here 
the percentage error of the exponent expansion with respect to the exact result 
(2.29) is plotted for various At, and compared with the approach of Rcf. [13]. The 
parameter choice, also taken from Ref. [13] corresponds to a sensible parameteriza- 
tion for interest rate markets. We adopt one year as unit of time, and we express 
the various parameters in this unit. The Euler approximation 

(x-x - p x (x )At) 2 " 



Pe(x, At\x ) = 



1 

2na 2 At 



exp 



2a 2 At 



(2.37) 
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is also reported for comparison. The inclusion of each successive order allows one 
to increase dramatically the accuracy of the approximation so that the third order 
expansion has basically a negligible error even for a sizable time step of order 6 
months. Remarkably, for the considered example, the third order expansion allows 
one to estimate the 10 years transition probability with a relative error of less than 
10 basis points (Fig. l-(b)). As illustrated in Fig. l-(a), in the present case the ap- 
proach of Ref. [13] provides a slightly poorer level of accuracy for n < 2. In addition, 
it generally produces more complicated mathematical expressions. Furthermore, as 
shown in Fig. l-(c), in the regime of small volatility the exponent expansion still 
provides accurate results while the performance of the approach of Ref. [13] de- 
grades. In fact, as anticipated, in the limit of small volatility (a x < 0.5%) the first 
order correction of the latter approach produces a negative transition probability 
signaling a break down of the scheme. 

2.2.2. Cox, Ingersoll and Ross diffusion 

The Vasicek model is probably too easy of a test case as the associated transition 
probability is Gaussian. In fact, since the exponent expansion has a leading term 
which is Gaussian, the higher powers in At just have to renormalize its average 
and variance in order to reproduce the exact result. It is interesting therefore to 
test the accuracy of the exponent expansion for a diffusion process that, while still 
integrable, generates a non-normal transition density. This is the case for the Feller's 
square root process [15] 



which is the basis of Cox, Ingersoll and Ross model for the instantaneous interest 
rate [3]. The exact transition probability is given by [3] 



where c = 2a/ [<jy ( 1 — cxp (— aAt))], q = 2ab/o-y — 1 > 0, u = cy cxp (— aAt), v = cy 
and I q is the modified Bessel function of the first kind of order q [16]. 

As explained in Sec. 2, since the volatility is not uniform, it is convenient to 
introduce the auxiliary process defined by Eq. (2.1), as X t = j(Y t ) = 2\fYt/a y . The 
X t process follows Eq. (2.2), with a x = 1 and 




(2.38) 




(2.39) 




x 



(2.40) 



and q = q + 1/2. In this case the effective potential reads 



t/ t \ 1 t \2 Q a 
V eS (x) = -fi x (x) -^2-4 



(2.41) 
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and the four terms of the exponent expansion in Eq. (2.8) are: 

„2 „2\ 



W (x,x Q ) = -glog h -(x 2 - x'q) 

Xq 4 



Wi(x,x ) 



2Ax 
a 2 



W 2 (x,x ) 
W 3 (x,x ) 



2Ax 2 
1 



Vx(xt) - Vx{xo) - q 2 

- -4) -aq(x- x Q ) 
[V eS (x) + V cS (x ) - 2W 1 (x, x 



2Ax 2 
3 



G(x t ) - G(x ) 



Ax 



W 2 {x,x ) + 



1 



[d x Ves{x) - d x V e s(x )] 



(2.42) 

(2.43) 
(2.44) 

(2.45) 



(2.46) 



Ax 2 > ^ 4Ax 3 

where Ax = x — x , d x V c ff(z) = q(l — q)/z 3 + a 2 z/A, and 

G(z) = U 2 z 5 \^ + 7 2 z + 2aPz -2^ + % 

and a = a 2 /8, = q(q — l)/2, 7 = — a(q + l)/2. Finally, going back to the original 
process, using Eq. (2.4), the transition probability reads: 

Py(y, At|j/ ) = Px{2y/y/<J Vl At\ 2y/yo~/(T y )/<Tyy/y . (2.47) 

The accuracy of the exponent expansion in this case is illustrated in Fig. 2. Simi- 
larly to the case of the Vasicek diffusion, the exponent expansion is characterized by 
a remarkably fast convergence by including successive terms of the approximation 
so that n — 3 provides already a virtually exact representation of the transition 
density, for At ~ lyrs. In this case, the approach of Ref. [13] performs slightly 
worse of the exponent expansion for n = 1, and slightly better for n = 2. However, 
also in this case the former breaks down for small values of the volatility, generating 
unphysical transition densities. 



2.2.3. Constant Elasticity of Variance diffusion 

As a last example we consider the Constant Elasticity of Variance model: 

dY t = a(b - Y t )dt + a y Y t p dW t (2.48) 

Here we consider for brevity only the case p > 1 and the transformation to a process 
with constant (unit) variance is X t = j(Y t ) = Y t ~ p /&y(p— 1) and gives, according 
to Eq. (2.3) 



Cl 



Hx{x) = — + c 2 x + C3XP-1 



(2.49) 



with ci =p/2(p-l), c 2 = a(p-l), and c 3 = -ab{p- l) p ^ p -^al /{p 1} . In this case 
the effective potential reads 



Cl 



C2 



Kff(z) = ~H X {X) - —r + - 



2x 2 2 2{p-l) 



^cax 1 ^-!) 



(2.50) 
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1 2 

Order of the Expansion (n) 



Fig. 2. Accuracy of the exponent expansion for the transition density of the Cox-Ingcrsoll-Ross 
model (2.38), for a = 0.0721, b = 0.219, a = 0.06665, and xq = 0.06. Bottom: comparison between 
the exponent expansion (stars) and the approximation of Ref. [13] (squares) for At = 0.5. The 
uniform error of the Euler approximation is also reported for comparison. Top: probability density 
function for At = 1.5 as a function of n: n = (dotted), n = 1 (long dashed), n = 2 (short dashed), 
n = 3 (continuous), Euler (dot-long dashed), exact (crosses). The inset is an enlargement of the 
region of the maximum. On this scale, the estimates for n = 2 and n = 3 still appear coincident. 



and the first three terms of the expansion are: 



W (x,x ) = cilog 



Wi(x,x ) 



yo 



a(p - 1) 



(2p-l)(y?-y 2 ) 



Vt 2(2p-l) 

1 y 2—1 ~ ^ \ 

26(^-1)^(7^ (xT^r -xj^) 
F(x) - F(x ) 



1 

2Ax 
1 



W 2 {x,x ) = [V*b{x) + V cfi (x ) - 2W!{x, x )] 



(2.51) 
(2.52) 
(2.53) 
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with 



z 3 3p - 



1) ^ 

-z p- 



1 



+ 2c 1 c 2 z + 



2cic 3 (p - 1) 



ZP- 



+ 



2C2C3(p - 1) -^+M,W. (2.54) 



p 3p- 2 

Similarly to the examples considered previously, also for the Constant Elasticity 
of Variance model we find a very fast convergence of the exponent expansion for 
At ~ 1, and a performance generally similar to the one of the approach of Ref. [13], 
for values of the volatility large enough. 



3. Pricing kernels for contingent claims 
3.1. Path integral formulation 

The exponent expansion can be generalized to obtain an approximation of the pric- 
ing kernels of 'standard' derivatives. This can be done by formulating the pricing 
problem within Feynman's path integral framework [6]. Here we indicate as 'stan- 
dard' any contingent claim written on the underlying, Y t , whose value at time t = 0, 
Vo, can be expressed as an expectation value of a functional of a certain type, namely 



V (At,y ) =E P(Y At )F[Y u ] y , 



where 



F[Y U ] = exp 



A/ 



duV F [Y u ] 



(3.1) 



(3.2) 



and P(Y/± t ) is a payout function. Above At is the time to expiry, and the expectation 
value is performed with respect to the probability measure defined by the diffusion 
for Y t that we assume of the form (1.1). European Vanilla options, zero coupon 
bonds, caplets, and floorets clearly belong to this family of contingent claims. In 
addition, other path-dependent derivatives, like barrier or Asian options can be 
expressed in this form (see e.g., Refs. [18,19]). 

Similarly to the case of the transition probability, it is in general convenient to 
introduce an auxiliary diffusion with constant volatility of the form (2.2) by means 
of the integral transformations (2.1). Then, the expectation value in (3.1) can be 
transformed in an integral over the distribution generated by such auxiliary diffusion 
by means of the usual Jacobian transformation (2.4). As a result, the value of the 
option can be in general written as: 



V (At,x ) = E P(X At )F[X u ] 



X 


-L 







dxP(x)K(x,At\x ) . 



(3.3) 



where D is the domain of the auxiliary process as defined by the relative stochastic 
differential equation (2.2), and K{x, At\xo) is the pricing kernel. The latter can be 
expressed in terms of a path integral as it follows [6] 

K(x, At\x ) = e- w °( x ' x °)$(x, At\x ) (3.4) 
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with 

px(At)=x I" /-At 

$>(x,At\xo) = / T>[x(u)]exp — / duL e s[x(u),x(u)] 

Jx(0)=xo JO 



(3.5) 



where Wq is given by Eq. (2.9) and the functional L e g [x(u), x(u)] is the effective 
Euclidean Lagrangian 

L cS [x(u),x(u)} = — ^i(w) 2 + V e s(x) (3.6) 

i0 ~x 

(x(u) = dx{u)/du) with the effective potential, V e g{x), defined as: 

V cS (x) = 7^2Vx{x) 2 + ^d xl i x {x) + V F (x) . (3.7) 

Finally, the measure V[x{u)\ is defined by discretizing each path x{u) connecting 
x(0) = xq and x(T) = x. This can be done by dividing the time interval [0,T] into 
P intervals so that x n — x{u n ) (u n = nT / P with n = 0, P), and by integrating 
the internal P — 1 variables x n over the domain D. The path integral J T>[x(u)] is 
then obtained as the limit for P — > oo of this procedure, namely 



/ 



p-i 

V[x{u)] = lim {2tto 2 x M)- p I 2 T\ / dx n . (3.8) 



It is well known form the physical sciences that the path integral $(x, At\ro) 
satisfies the partial differential equation Eq. (2.25) [8]. Note that this is consistent 
with the fact that, by definition, for Vf(x) = the pricing kernel coincides with 
the transition density of the underlying diffusion process for X t . In particular, as 
observed in Sec. 2.1, one can use Eq. (2.25) to derive the exponent expansion for 
&(x, Ax\ro) using the trial form (2.27). As a result, the same expressions Eqs.(2.16- 
2.18) derived for the transition density hold true also for the pricing kernel, provided 
that the effective potential (3.7) replaces the one in Eq. (2.19). 



3.1.1. Correspondence with Quantum Mechanics 

It is interesting to note that the path integral formulation of the pricing kernel (3.5) 
is mathematically equivalent to the quantum mechanical description of the thermo- 
dynamic properties of an ideal gas of particles moving in the potential hV c s{x) (h is 
the reduced Planck's constant giving the correct energy dimensions). The complete 
correspondence is obtained by identifying a 2 — ► h/m and At — > h/ksT where m is 
the mass of the particle, T is the temperature, and ks is the Boltzmann constant. 
With this prescription, &(x, At\x ) becomes the so-called single particle density ma- 
trix, and the results of Makri and Miller [5] can be readily recovered. In addition, 
it is straightforward to show using Eqs. (2.23) that the exponent expansion of its 
diagonal elements, $>(xo,At\xo), are consistent with the so called Wigncr expan- 
sion for the high-temperature limit. Finally, performing the analytic continuation 
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x 



Fig. 3. Accuracy of the exponent expansion for the pricing kernel Eq. (3.4) of the Vasicck model 
(2.28), for a = 0.0717, b = 0.261, and xq = 0.1. Symbols as in Fig. 2. The inset shows the maximum 
absolute error as a function of the order of the expansion n. 

known as Wick rotation At — ► iht allows one to obtain the single-particle quantum 
propagator. In this case (2.25) becomes the celebrated Schrodinger equation. 

This correspondence provides an alternative justification of the choice of the 
exponential ansatz in Eq. (2.5). Indeed, this is the form in which can be expressed 
in general the quantum mechanical propagator or the single particle density matrix 
[8]. Furthermore, it has been shown for the quantum problem [17] that the exponent 
expansion up to third order in At and first order in h/m can be derived starting 
from the short time semiclassical propagator obtained [20] through a saddle point 
analysis of the limit h/m — ► 0. Indeed, it can be shown that the second order 
correction W% (x,Xo), Eq. (2. 17) , is the so-called van-Vleck determinant of the saddle 
point expansion. This explains why the accuracy of the present scheme is preserved 
in the corresponding regime of low volatility, as it was anticipated in Sec. 2, and 
illustrated in Sec. 2. 2. 



3.1.2. An example: Zero Coupon Bond 

In this Section we illustrate the prescriptions outlined above by applying the ex- 
ponent expansion to the calculation of a zero coupon bond within the Vasicek and 
Cox-Ingersoll-Ross models. The zero coupon bond is a financial instrument that 
provides at time t — At a payout of one unit of a certain notional. Its value at time 
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P(0,Ai) = E 



exp - 



A/ 



duX-, 



(3.9) 



which is of the standard form given by Eqs. (3.1) and (3.2), with = X u , 

and P[X At ] = 1. 

As a result, the exponent expansion for the kernel Eq. (3.4) can be easily derived 
giving for the Vasicek model 

x + x 



W 1 (x,x ) = W?(x,x ) + 
W 2 (x 1 x ) =W°(x,x ) 

W 3 (x,x ) = W§(x,x ) - — 



(3.10) 
(3.11) 
(3.12) 



where Wf(x,xo) are the expressions obtained for the transition probability 
Eqs. (2.33)-(2.36) of Sec. 2.2.1. For the Cox-Ingersoll-Ross model instead we get: 

(3.13) 

(3.14) 



Wi(x,x ) = W°(x,x ) + — (x 2 +xl + xx ) 



W 2 (x,x )=W°(x,x )- — 

with Wi(x,x ) given by Eqs. (2.42)-(2.44), and W§(x,xo) related to the previous 
quantities as in Eq. (2.45) with 



(3.15) 



G°(z) as in Eq. (2.46), and d z V e g(z) — d z V® s (z) + <J 2 z/2. The exponent expansion 
for the pricing kernel can be compared with the exact results that can be shown to 
read for the Vasicek model (2.28) 



K cx (x, At\,x ) 



exp [(x - x )/a - At(b - a 2 /2a 2 )] 
(2tt ( t 2 ) 1 /2 



exp 



[(x - b + o 2 ja 2 ) e- aAt -{y-b + a 2 /a 2 ^ 



with a given by Eq. (2.30), and 



K e x(x, At\,xo) = - exp 
x 



2a 2 



(x 2 -x 2 ) + (2ab- ^)log ; 



(3.16) 



' 4 



"fy/xx a e 



a 2 bAt/a 2 



exp 



7 
4ct 2 



(x 2 + x 2 ) coth [7 At/2] 



In 



XX "f 



2a 2 sinh [7 At/2] 



2a 2 sinh [7 At/2] 

(3.17) 



with 7 = Va 2 + 2cr 2 , for the Cox-Ingersoll-Ross one. As illustrated in Fig. 3, sim- 
ilarly to the case of the transition probability, the exponent expansion provides a 
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Fig. 4. Zero Coupon Bond for the Cox-Ingersoll-Ross model. Parameters and symbols as in Fig. 2. 

remarkably good, and fast converging approximation of the exact pricing kernel for 
financially sensible parameterizations, and for a sizable value of the time step At. 

Finally, the zero coupon bond can be obtained by numerical integration of the 
pricing kernel according to Eq. (3.3). The corresponding results are shown in Fig. 4 
confirming once more the quality of the approximation. In a similar fashion, one 
can obtain systematic approximations for caplets, floorets, and other simple interest 
rate derivatives whose value depends only on the value of the instantaneous short 
rate at time t = At. It is also possible to generalize this approach to path dependent 
contingent claims, like Asian options. 



0.5 I i i i i | i i i i 




4. Extending the time-step: path integral Monte Carlo methods 

For an extended time interval T, the calculation of the transition density or, more in 
general, of the pricing kernel (3.4) can be performed by discretizing the path integral 
(3.5), and taking the limit of large number of time slices (P — > oo) according to the 
standard Trotter product formula: 

p-i 



K(x,At\x ) ~ {2ira 2 x At)- p / 2 e- w °( x ' X0 *> TT f dx 

n=l jD 

k=P 



x exp 



2alAt 



k=P 



(4.1) 



Xk-lf - — ^2 ( V efl( x k) + Kff(Zfc-l)) 
fe=l fe=l 

with At = T/P, xp = x, and the effective potential given by Eq. (3.7). It is 
worth remarking that, for the case of the transition density, by interpreting the 
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latter equation as the Chapman-Kolmogorov property of Markov processes [7] one 
obtains the following approximation of the short-time propagator 



-^Trotter At\x ) 



,-W h 



exp 



(x - Xq) 2 

2alAt 



— \V eS {x) + V cS (x ) 



(4.2) 

However, in contrast to the n = 1 exponent expansion, the latter expression is not 
strictly correct up to order At, and only in the limit P — > oo the difference becomes 
negligible. 

In general, to obtain an accurate result for the pricing kernel for an extended time 
period T one has to increase the number of time slices, or Trotter number P, until 
convergence is achieved. By replacing the Trotter formula with the improved short- 
time kernel obtained through the exponent expansion (2.5) one achieves a faster 
convergence with the Trotter number, thus significantly reducing the computational 
burden. In this case the finite-time expression of the pricing kernel reads 



p-i 

K(x,T\x )c(2na 2 x At)- p / 2 ]J / 

~, — i J L 



x exp 



1 P P 

2X7 ^2( X k - Xk-lf - ^2 W ( X k,Xk-l,&t) 



2olAt 



k=l 



fc=l 



(4.3) 



with W(x, x , At) given by Eq. (2.8). 

Equation (4.3) allows one to obtain the transition density or the pricing kernel 
for a standard derivative through the calculation of a multidimensional integral over 
the variables x\, . . . ,xp-\. The latter integration is ideally suited for Monte Carlo 
methods either in real, or in Fourier space [21], the specific choice depending on 
the particular problem at hand. In addition, importance-sampling schemes, e.g., by 
means of the Metropolis algorithm [22] , can be easily applied in order to reduce the 
computation time. However, in order not to introduce a systematic bias in the result 
a particular attention has to be paid in order to sample accurately the configuration 
space. 

The most straightforward way to perform a Monte Carlo quadrature of Eq. (4.3) , 
is to realize that a simple Markov chain x = (x\, . . . , Xp-i) 



Xn-l + (T x \fAt 



(4.4) 



with Z n , sampled from a standard normal distribution, generates an ensemble of 
walkers distributed according to 



p(xi, . . .,xp-i\x ) = {2na- 2 x At) 



-(P-l)/2 



exp 



I 



k=l 



(4.5) 



As a result, the pricing kernel (4.3) can be obtained as the average over the random 
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paths generated according to Eq. (4.4) of the following estimator: 



0{x lXp =x) = {2iralAty 1/2 exp 



! p 

— (X - Ip-l) 2 - ^ W ( X k, Xk-l) 



lalAt 



k=l 



(4.6) 

A remarkable property of the path integral approach is that K(x,At\xo) for any 
final point x can be evaluated with a single Monte Carlo simulation by appropri- 
ately reweighting the accumulated estimator. In fact the distribution of walkers 
p(xi, . . . , xp-i\xo) is independent of the final point xp so that K(x',At\xo) can 
be calculated by averaging 0(x, xp = x'). In addition, the latter quantity can be 
efficiently obtained by means of the following reweighting procedure 

W(x',x) 



0(x, xp = x') = 0(x, xp = x) 



W(ar,x) 



with 



W(x, x) = exp 



1 



2cr?At 



(x — xp-i) 2 — W(x, Xp_\) 



(4.7) 



(4.8) 



Expectation values of the form (3.3) on a time horizon T can be calculated by 
integrating over the final variable giving: 

p 

V {T,x )= [ dxP(x)K(x,T\x a )~(27ra 2 x At)- p / 2 t[ [ dx n P(x P ) 
Jd u=1 Jd 

1 P P 



exp 



2alAt 



k=i 



fe=i 



(4.9) 



This can be simulated by extending the Markov chain (4.4) up to step and 
accumulating the estimator 



P(x,xp = x 1 ) = P(x) exp 



^IF(x fe ,a: fe _i) 



k=i 



(4.10) 



Remarkably, within the path integral approach, the sensitivities of such expectation 
values with respect to the model parameters (the so-called Greeks) can be computed 
in the same Monte Carlo simulation, thus avoiding any numerical differentiation. 
Indeed, indicating with a generic parameter, under quite general conditions [23], 
one has 

d e V (T,x ,6) = [ dx[Kg(x,T\x a )dgP(x,e)+P(x,9)deKe(x,T\x a )} . (4.11) 
Jd 

As a result the sensitivity dgVo(T, xo, 0) can be calculated by means of the estimator: 

p 

^W{x k ,x k ^ 



Q(x, xp = x') = exp 



fe=i 



(deP + PdglogKg) . 



(4.12) 



Higher order sensitivities can be obtained in a similar fashion. 
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Fig. 5. Convergence with the Trotter number P of the normalization, and of the first 4 moments 
of the 40 years transition density for the CIR model. Parameters as in Fig. 2. In the first top panels 
the relative error with respect to the exact result is reported. Lines are quadratic fits. 



The convergence with the Trotter number P of the path integral Monte Carlo 
estimates is illustrated in Fig. 5 for the calculation of the first five moments of the 
T = 40yrs transition probability of the Cox-Ingersoll-Ross model (2.38). The finite 
P estimates converge very rapidly with 1/P. In particular, for the case considered, 
P = 20 already provides estimates in agreement with the exact result within statis- 
tical uncertainties. In general, as also shown in the figure, a convenient indicator of 
the convergence is the zeroth-moment or normalization of the distribution. The cal- 
culation of this quantity allows in general to assess the level of convergence without 
performing a complete scaling with P, thus saving computational time. 



5. Conclusions 

The exponent expansion is an approximation of the quantum mechanical propagator 
known in physical chemistry [5,17]. We have generalized this approach to produce 
closed-form approximation of the transition probability of arbitrary non linear pro- 
cesses, and we have shown that it produces very accurate results for integrable 
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diffusions of financial interest, like the Vasicek and the Cox-Ingersoll-Ross models. 
In contrast to previously developed approximation schemes that share a similar 
rationale [13], the exponent expansion always generates positive definite transition 
probabilities, and remains stable also in the limit of low volatility. 

By introducing a path integral framework we have generalized the exponent 
expansion to the calculation of pricing kernels of financial derivatives, and we have 
shown how to obtain accurate approximations for the price of simple contingent 
claims. We have also shown how the exponent expansion can be naturally used to 
increase the efficiency of path integral Monte Carlo simulations. The latter allow 
one to extend the calculations to arbitrarily large time steps, and to efficiently 
calculate the Greeks of contingent claims, avoiding any numerical differentiation. A 
systematic study of the efficiency of this approach for the pricing of exotic derivatives 
will be the object of future investigations. 

The exponent expansion can be generalized to multifactor diffusion processes, 
and to time dependent drift and volatility functions. This would allow one to in- 
crease the efficiency of numerical simulations of a larger family of diffusion processes 
relevant for Financial applications, including local volatility or Libor Market Mod- 
els. Work is currently in progress along this line of research. 
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